BioGSP Tutorial: Spatial Analysis with Toy CODEX Data

Introduction

This tutorial demonstrates the BioGSP (Biological Graph Signal Processing) package for analyzing spatial patterns in multiplexed imaging data. We’ll use the built-in toy CODEX dataset to showcase the package’s capabilities including:

• Multiple kernel family options: Mexican hat, Meyer, and heat kernel families
• Spatial pattern decomposition: Multi-scale analysis of cell type distributions
• Graph construction and analysis: Eigenvalue analysis and graph properties
• Cross-cell type correlation: Advanced correlation measures
• Comprehensive visualization: Filter banks, decompositions, and spatial patterns

Load Required Libraries

# Load required libraries
library(BioGSP)
library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(knitr)

# Verify BioGSP package version
cat("BioGSP Package version:", as.character(packageVersion("BioGSP")), "\n")

## BioGSP Package version: 1.0.0

BioGSP_Tutorial.R

Load and Explore the Toy Dataset

# Load the built-in toy CODEX dataset
data(codex_toy_data)

# Examine the dataset structure
str(codex_toy_data)

## 'data.frame':    18604 obs. of  5 variables:
##  $ cellLabel  : chr  "ROI_0_BCL6BCell_1" "ROI_0_BCL6BCell_2" "ROI_0_BCL6BCell_3" "ROI_0_BCL6BCell_4" ...
##  $ Y_cent     : num  64.5 71.2 78.3 50.8 52.6 ...
##  $ X_cent     : num  112 82.6 91.4 92.6 98.8 ...
##  $ Annotation5: chr  "BCL6- B Cell" "BCL6- B Cell" "BCL6- B Cell" "BCL6- B Cell" ...
##  $ ROI_num    : chr  "ROI_0" "ROI_0" "ROI_0" "ROI_0" ...

head(codex_toy_data)

##                 cellLabel   Y_cent    X_cent  Annotation5 ROI_num
## ROI_0.1 ROI_0_BCL6BCell_1 64.52399 111.96648 BCL6- B Cell   ROI_0
## ROI_0.2 ROI_0_BCL6BCell_2 71.22437  82.56243 BCL6- B Cell   ROI_0
## ROI_0.3 ROI_0_BCL6BCell_3 78.33008  91.44301 BCL6- B Cell   ROI_0
## ROI_0.4 ROI_0_BCL6BCell_4 50.83888  92.60675 BCL6- B Cell   ROI_0
## ROI_0.5 ROI_0_BCL6BCell_5 52.64693  98.76092 BCL6- B Cell   ROI_0
## ROI_0.6 ROI_0_BCL6BCell_6 56.20720  91.39930 BCL6- B Cell   ROI_0

# Summary of cell types
cell_summary <- codex_toy_data %>%
  group_by(Annotation5) %>%
  summarise(
    Count = n(),
    Mean_X = round(mean(X_cent), 2),
    Mean_Y = round(mean(Y_cent), 2),
    SD_X = round(sd(X_cent), 2),
    SD_Y = round(sd(Y_cent), 2)
  ) %>%
  arrange(desc(Count))

kable(cell_summary, caption = "Summary of Cell Types in Toy CODEX Data")

Summary of Cell Types in Toy CODEX Data

Annotation5	Count	Mean_X	Mean_Y	SD_X	SD_Y
CD4 T	4092	45.64	45.52	21.49	23.61
BCL6- B Cell	3719	55.43	45.64	22.96	20.28
CD8 T	3346	47.38	48.53	23.91	27.34
DC	2233	50.17	46.15	23.18	24.00
CD4 Treg	1490	53.46	40.06	23.25	19.53
M1	1490	46.59	54.36	21.19	22.62
BCL6+ B Cell	931	36.53	51.97	21.03	23.56
Endothelial	746	51.29	45.41	22.90	19.33
M2	370	56.67	37.24	28.21	19.59
Myeloid	186	47.60	39.69	21.09	20.63
Other	1	107.12	46.12	NA	NA

# Focus on BCL6- B Cell and CD4 T for detailed analysis
target_cells <- c("BCL6- B Cell", "CD4 T")
cat("\nTarget cell types for SGWT analysis:\n")

## 
## Target cell types for SGWT analysis:

cat("BCL6- B Cell:", sum(codex_toy_data$Annotation5 == "BCL6- B Cell"), "\n")

## BCL6- B Cell: 3719

cat("CD4 T cells:", sum(codex_toy_data$Annotation5 == "CD4 T"), "\n")

## CD4 T cells: 4092

# ROI summary
cat("\nROI summary:\n")

## 
## ROI summary:

roi_summary <- table(codex_toy_data$ROI_num)
print(roi_summary)

## 
##  ROI_0  ROI_1 ROI_10 ROI_11 ROI_12 ROI_13 ROI_14 ROI_15  ROI_2  ROI_3  ROI_4 
##    952    945   1315   1045   1190   1174   1497   1316   1155   1422   1097 
##  ROI_5  ROI_6  ROI_7  ROI_8  ROI_9 
##   1420   1059    958   1114    945

cat("\nCell types by ROI:\n")

## 
## Cell types by ROI:

cross_tab <- table(codex_toy_data$ROI_num, codex_toy_data$Annotation5)
print(cross_tab[, c("BCL6- B Cell", "CD4 T", "CD8 T", "DC", "M1")])

##         
##          BCL6- B Cell CD4 T CD8 T  DC  M1
##   ROI_0           190   209   171 114  76
##   ROI_1           189   208   170 113  76
##   ROI_10          263   289   237 158 105
##   ROI_11          209   230   188 125  84
##   ROI_12          238   262   214 143  95
##   ROI_13          235   258   211 141  94
##   ROI_14          299   329   269 180 120
##   ROI_15          263   290   237 158 105
##   ROI_2           231   254   208 139  92
##   ROI_3           284   313   256 171 114
##   ROI_4           219   241   197 132  88
##   ROI_5           284   312   256 170 114
##   ROI_6           212   233   190 127  85
##   ROI_7           191   211   172 115  77
##   ROI_8           223   245   200 134  89
##   ROI_9           189   208   170 113  76

BioGSP_Tutorial.R

Visualize Spatial Distribution

# Color mapping for visualization
fill_colors <- c("#e07329", "#16964a", "#2958a8", "#dd2246", "#703594", 
                 "#d0bd2a", "#69bd48", "#2297b1", "#b9e4f3", "#ae5c71", "#ab2170")
color_mapping <- data.frame(
  Annotation5 = c("BCL6- B Cell", "BCL6+ B Cell", "CD4 T", "CD4 Treg", "CD8 T", "DC", "Endothelial", "M1", "M2", "Myeloid", "Other"),
  Color = fill_colors
)
color_mapping_vector <- setNames(color_mapping$Color, color_mapping$Annotation5)

# Plot overall spatial distribution
p_overview <- ggplot(codex_toy_data, aes(x = X_cent, y = Y_cent, color = Annotation5)) +
  geom_point(size = 0.8, alpha = 0.7) +
  scale_color_manual(values = color_mapping_vector) +
  labs(title = "Toy CODEX Spatial Cell Distribution",
       x = "X Coordinate", y = "Y Coordinate") +
  theme_minimal() +
  scale_y_reverse() +
  guides(color = guide_legend(override.aes = list(size = 2)))

print(p_overview)

# Show distribution by ROI
p_by_roi <- ggplot(codex_toy_data, aes(x = X_cent, y = Y_cent, color = Annotation5)) +
  geom_point(size = 0.8, alpha = 0.7) +
  facet_wrap(~ROI_num, scales = "free") +
  scale_color_manual(values = color_mapping_vector) +
  labs(title = "Spatial Distribution by ROI",
       x = "X Coordinate", y = "Y Coordinate") +
  theme_minimal() +
  scale_y_reverse() +
  theme(legend.position = "bottom")

print(p_by_roi)

# Highlight target cell types
target_data <- codex_toy_data %>% filter(Annotation5 %in% target_cells)
p_target <- ggplot(codex_toy_data, aes(x = X_cent, y = Y_cent)) +
  geom_point(color = "lightgray", size = 0.3, alpha = 0.3) +
  geom_point(data = target_data, aes(color = Annotation5), size = 1.2) +
  facet_wrap(~ROI_num, scales = "free") +
  scale_color_manual(values = color_mapping_vector[target_cells]) +
  labs(title = "BCL6- B Cell and CD4 T Cell Distribution by ROI",
       x = "X Coordinate", y = "Y Coordinate") +
  theme_minimal() +
  scale_y_reverse()

print(p_target)

BioGSP_Tutorial.R

Kernel Family Comparison

cat("Demonstrating different kernel family options available in SGWT...\n")

## Demonstrating different kernel family options available in SGWT...

# Compare different kernel families
kernel_comparison <- compare_kernel_families(x_range = c(0, 3), 
                                            scale_param = 1, 
                                            plot_results = TRUE)

# Show available kernel family types
kernel_types <- c("mexican_hat", "meyer", "heat")
cat("\nAvailable kernel family types:\n")

## 
## Available kernel family types:

for (i in seq_along(kernel_types)) {
  cat(i, ". ", kernel_types[i], "\n", sep = "")
}

## 1. mexican_hat
## 2. meyer
## 3. heat

# Show what each kernel family provides
cat("\nEach kernel family provides:\n")

## 
## Each kernel family provides:

cat("- mexican_hat: Gaussian scaling + LoG-style wavelet\n")

## - mexican_hat: Gaussian scaling + LoG-style wavelet

cat("- meyer: Smooth cosine low-pass + band-pass transitions\n") 

## - meyer: Smooth cosine low-pass + band-pass transitions

cat("- heat: Exponential decay scaling + derivative-like wavelet\n")

## - heat: Exponential decay scaling + derivative-like wavelet

# Display first few rows of comparison
cat("\nKernel values at different eigenvalues:\n")

## 
## Kernel values at different eigenvalues:

print(head(kernel_comparison, 8))

##            x mexican_hat_scaling mexican_hat_wavelet meyer_scaling
## 1 0.00000000           1.0000000        0.0000000000             1
## 2 0.01507538           0.9998864        0.0002272412             1
## 3 0.03015075           0.9995456        0.0009086548             1
## 4 0.04522613           0.9989778        0.0020433121             1
## 5 0.06030151           0.9981835        0.0036296666             1
## 6 0.07537688           0.9971632        0.0056655569             1
## 7 0.09045226           0.9959176        0.0081482106             1
## 8 0.10552764           0.9944474        0.0110742486             1
##   meyer_wavelet heat_scaling heat_wavelet
## 1             0    1.0000000   0.00000000
## 2             0    0.9850377   0.01484981
## 3             0    0.9702992   0.02925525
## 4             0    0.9557813   0.04322629
## 5             0    0.9414806   0.05677270
## 6             0    0.9273939   0.06990406
## 7             0    0.9135179   0.08262976
## 8             0    0.8998496   0.09495900

BioGSP_Tutorial.R

Create Binned Data for SGWT Analysis

# Function to create binned data for SGWT analysis
create_binned_data <- function(cell_data, cell_type, x_bins = 30, y_bins = 30) {
  # Filter for specific cell type
  filtered_data <- cell_data %>% filter(Annotation5 == cell_type)
  
  if (nrow(filtered_data) == 0) {
    stop(paste("No cells found for cell type:", cell_type))
  }
  
  # Create bins
  filtered_data <- filtered_data %>%
    mutate(
      x_bin = cut(X_cent, breaks = x_bins, labels = FALSE),
      y_bin = cut(Y_cent, breaks = y_bins, labels = FALSE)
    )
  
  # Count cells in each bin
  bin_counts <- filtered_data %>%
    group_by(x_bin, y_bin) %>%
    summarise(cell_count = n(), .groups = 'drop')
  
  # Create complete grid
  complete_grid <- expand.grid(x_bin = 1:x_bins, y_bin = 1:y_bins)
  
  # Merge with counts (fill missing with 0)
  binned_data <- complete_grid %>%
    left_join(bin_counts, by = c("x_bin", "y_bin")) %>%
    mutate(cell_count = ifelse(is.na(cell_count), 0, cell_count))
  
  # Calculate proportions
  binned_data$cell_proportion <- binned_data$cell_count / max(binned_data$cell_count)
  
  return(binned_data)
}

# Create binned data for both cell types in ROI_1
# Filter for ROI_1 for focused analysis
roi1_data <- codex_toy_data %>% filter(ROI_num == "ROI_1")

bcl6nb_binned <- create_binned_data(roi1_data, "BCL6- B Cell")
cd4t_binned <- create_binned_data(roi1_data, "CD4 T")

cat("Binned data created successfully!\n")

## Binned data created successfully!

cat("BCL6- B Cell bins with cells:", sum(bcl6nb_binned$cell_count > 0), "out of 900\n")

## BCL6- B Cell bins with cells: 157 out of 900

cat("CD4 T bins with cells:", sum(cd4t_binned$cell_count > 0), "out of 900\n")

## CD4 T bins with cells: 171 out of 900

# Visualize binned data
plot_binned_data <- function(binned_data, cell_type, color_val) {
  ggplot(binned_data, aes(x = x_bin, y = y_bin, fill = cell_proportion)) +
    geom_tile() +
    scale_fill_gradient2(low = "white", high = color_val, name = "Proportion") +
    labs(title = paste(cell_type, "Cell Distribution (30x30 bins)"),
         x = "X Bin", y = "Y Bin") +
    theme_minimal() +
    scale_y_reverse() +
    coord_equal()
}

p_bcl6nb_bins <- plot_binned_data(bcl6nb_binned, "BCL6- B Cell", color_mapping_vector["BCL6- B Cell"])
p_cd4t_bins <- plot_binned_data(cd4t_binned, "CD4 T", color_mapping_vector["CD4 T"])

# Show side by side
bins_combined <- p_bcl6nb_bins | p_cd4t_bins
print(bins_combined)

BioGSP_Tutorial.R

SGWT Analysis

Prepare Data for SGWT

# Prepare data for SGWT analysis
prepare_sgwt_data <- function(binned_data) {
  sgwt_data <- binned_data %>%
    mutate(
      x = x_bin,
      y = y_bin,
      signal = cell_proportion
    ) %>%
    select(x, y, signal)
  
  return(sgwt_data)
}

# Prepare data for both cell types
bcl6nb_sgwt_data <- prepare_sgwt_data(bcl6nb_binned)
cd4t_sgwt_data <- prepare_sgwt_data(cd4t_binned)

cat("SGWT data prepared:\n")

## SGWT data prepared:

cat("BCL6nB data points:", nrow(bcl6nb_sgwt_data), "\n")

## BCL6nB data points: 900

cat("CD4T data points:", nrow(cd4t_sgwt_data), "\n")

## CD4T data points: 900

cat("Signal range BCL6nB:", round(range(bcl6nb_sgwt_data$signal), 4), "\n")

## Signal range BCL6nB: 0 1

cat("Signal range CD4T:", round(range(cd4t_sgwt_data$signal), 4), "\n")

## Signal range CD4T: 0 1

BioGSP_Tutorial.R

Apply SGWT Analysis

cat("Applying SGWT to BCL6nB cell distribution...\n")

## Applying SGWT to BCL6nB cell distribution...

# Apply SGWT to BCL6nB data
sgwt_bcl6nb <- SGWT(data.in = bcl6nb_sgwt_data,
                    signal = "signal",
                    k = 15,  # Number of nearest neighbors
                    J = 3,   # Number of scales
                    scaling_factor = 2,
                    kernel_type = "mexican_hat",  # Unified kernel family
                    laplacian_type = "normalized",
                    k_fold = 2,
                    return_all = TRUE)

## Building graph from spatial coordinates...
## Computing Laplacian and eigendecomposition...
## Auto-generated scales: 0.4886, 0.2443, 0.1221 
## Performing SGWT decomposition...
## Reconstruction RMSE: 0.11396

cat("SGWT analysis completed for BCL6nB\n")

## SGWT analysis completed for BCL6nB

cat("Reconstruction RMSE:", round(sgwt_bcl6nb$reconstruction_error, 6), "\n")

## Reconstruction RMSE: 0.11396

# Apply SGWT to CD4T data
cat("\nApplying SGWT to CD4T cell distribution...\n")

## 
## Applying SGWT to CD4T cell distribution...

sgwt_cd4t <- SGWT(data.in = cd4t_sgwt_data,
                  signal = "signal",
                  k = 15,
                  J = 3,
                  scaling_factor = 2,
                  kernel_type = "mexican_hat",
                  laplacian_type = "normalized",
                  k_fold = 2,
                  return_all = TRUE)

## Building graph from spatial coordinates...
## Computing Laplacian and eigendecomposition...
## Auto-generated scales: 0.4886, 0.2443, 0.1221 
## Performing SGWT decomposition...
## Reconstruction RMSE: 0.149172

cat("SGWT analysis completed for CD4T\n")

## SGWT analysis completed for CD4T

cat("Reconstruction RMSE:", round(sgwt_cd4t$reconstruction_error, 6), "\n")

## Reconstruction RMSE: 0.149172

# Energy analysis
energy_bcl6nb <- sgwt_energy_analysis(sgwt_bcl6nb)
energy_cd4t <- sgwt_energy_analysis(sgwt_cd4t)

# Combine energy data for comparison
energy_bcl6nb$cell_type <- "BCL6nB"
energy_cd4t$cell_type <- "CD4T"
energy_combined <- rbind(energy_bcl6nb, energy_cd4t)

# Plot energy distribution
p_energy <- ggplot(energy_combined, aes(x = scale, y = energy_ratio, fill = cell_type)) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_fill_manual(values = c("BCL6nB" = color_mapping_vector["BCL6nB"],
                               "CD4T" = color_mapping_vector["CD4T"])) +
  labs(title = "Energy Distribution Across SGWT Scales",
       x = "Scale", y = "Energy Ratio", fill = "Cell Type") +
  theme_minimal()

print(p_energy)

BioGSP_Tutorial.R

Compare Different Kernel Families

cat("Testing different kernel families on BCL6nB data...\n")

## Testing different kernel families on BCL6nB data...

# Extract components for testing
eigenvalues <- sgwt_bcl6nb$graph_info$eigenvalues
eigenvectors <- sgwt_bcl6nb$graph_info$eigenvectors
signal_vector <- bcl6nb_sgwt_data$signal

# Test different kernel families
kernel_types_test <- c("mexican_hat", "meyer", "heat")
kernel_results <- list()

for (kernel_type in kernel_types_test) {
  cat("Testing", kernel_type, "kernel family...\n")
  
  result <- sgwt_forward(signal = signal_vector,
                        eigenvectors = eigenvectors,
                        eigenvalues = eigenvalues,
                        scales = sgwt_bcl6nb$parameters$scales,
                        kernel_type = kernel_type)
  
  reconstructed <- sgwt_inverse(result)
  rmse <- sqrt(mean((signal_vector - reconstructed)^2))
  
  kernel_results[[kernel_type]] <- list(
    result = result,
    reconstructed = reconstructed,
    rmse = rmse
  )
  
  cat("RMSE with", kernel_type, "kernel:", round(rmse, 6), "\n")
}

## Testing mexican_hat kernel family...
## RMSE with mexican_hat kernel: 0.11396 
## Testing meyer kernel family...
## RMSE with meyer kernel: 0.11527 
## Testing heat kernel family...
## RMSE with heat kernel: 0.113988

# Create comparison table
kernel_comparison_table <- data.frame(
  Kernel_Family = names(kernel_results),
  RMSE = sapply(kernel_results, function(x) round(x$rmse, 6)),
  Description = c("Gaussian scaling + LoG wavelet", 
                  "Smooth cosine transitions", 
                  "Exponential decay + derivative"),
  stringsAsFactors = FALSE
)

kable(kernel_comparison_table, 
      caption = "Reconstruction Error Comparison Across Different Kernel Families")

Reconstruction Error Comparison Across Different Kernel Families

	Kernel_Family	RMSE	Description
mexican_hat	mexican_hat	0.113960	Gaussian scaling + LoG wavelet
meyer	meyer	0.115270	Smooth cosine transitions
heat	heat	0.113988	Exponential decay + derivative

BioGSP_Tutorial.R

SGWT Decomposition Visualization

# Function to visualize SGWT components
plot_sgwt_components <- function(sgwt_result, data_in, cell_type, color_val) {
  coefficients <- sgwt_result$decomposition$coefficients
  plot_data <- data_in
  
  plots <- list()
  
  # Original signal
  plot_data$original <- sgwt_result$original_signal
  p_orig <- ggplot(plot_data, aes(x = x, y = y, fill = original)) +
    geom_tile() +
    scale_fill_gradient2(low = "white", high = color_val) +
    labs(title = paste(cell_type, "- Original")) +
    theme_void() + scale_y_reverse() + coord_equal() +
    theme(legend.position = "none")
  plots[["original"]] <- p_orig
  
  # Scaling function
  plot_data$scaling <- Re(as.vector(coefficients[[1]]))
  p_scaling <- ggplot(plot_data, aes(x = x, y = y, fill = scaling)) +
    geom_tile() +
    scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
    labs(title = "Scaling (Low-freq)") +
    theme_void() + scale_y_reverse() + coord_equal() +
    theme(legend.position = "none")
  plots[["scaling"]] <- p_scaling
  
  # Wavelet coefficients
  for (i in 1:min(3, length(coefficients) - 1)) {
    coeff_name <- paste0("wavelet_", i)
    plot_data[[coeff_name]] <- Re(as.vector(coefficients[[i + 1]]))
    
    p_wavelet <- ggplot(plot_data, aes_string(x = "x", y = "y", fill = coeff_name)) +
      geom_tile() +
      scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
      labs(title = paste("Wavelet", i)) +
      theme_void() + scale_y_reverse() + coord_equal() +
      theme(legend.position = "none")
    
    plots[[coeff_name]] <- p_wavelet
  }
  
  # Reconstructed signal
  plot_data$reconstructed <- sgwt_result$reconstructed_signal
  p_recon <- ggplot(plot_data, aes(x = x, y = y, fill = reconstructed)) +
    geom_tile() +
    scale_fill_gradient2(low = "white", high = color_val) +
    labs(title = "Reconstructed") +
    theme_void() + scale_y_reverse() + coord_equal() +
    theme(legend.position = "none")
  plots[["reconstructed"]] <- p_recon
  
  return(plots)
}

# Plot BCL6nB decomposition
bcl6nb_plots <- plot_sgwt_components(sgwt_bcl6nb, bcl6nb_sgwt_data, "BCL6nB", 
                                     color_mapping_vector["BCL6nB"])
bcl6nb_combined <- wrap_plots(bcl6nb_plots, ncol = 3)
print(bcl6nb_combined)

# Plot CD4T decomposition
cd4t_plots <- plot_sgwt_components(sgwt_cd4t, cd4t_sgwt_data, "CD4T",
                                   color_mapping_vector["CD4T"])
cd4t_combined <- wrap_plots(cd4t_plots, ncol = 3)
print(cd4t_combined)

BioGSP_Tutorial.R

Cross-Cell Type Correlation Analysis

cat("Performing cross-correlation analysis between BCL6nB and CD4T...\n")

## Performing cross-correlation analysis between BCL6nB and CD4T...

# Calculate Graph Cross-Correlation (GCC)
gcc_result <- Cal_GCC(
  data.in = data.frame(
    BCL6nB = bcl6nb_sgwt_data$signal,
    CD4T = cd4t_sgwt_data$signal
  ),
  knee = c(length(eigenvalues) %/% 4),
  signal1 = "BCL6nB",
  signal2 = "CD4T",
  eigenvector = eigenvectors
)

# Cosine similarity and Pearson correlation for comparison
cosine_sim <- cosine_similarity(bcl6nb_sgwt_data$signal, cd4t_sgwt_data$signal)
pearson_cor <- cor(bcl6nb_sgwt_data$signal, cd4t_sgwt_data$signal)

cat("Graph Cross-Correlation between BCL6nB and CD4T:", round(gcc_result, 4), "\n")

## Graph Cross-Correlation between BCL6nB and CD4T: 0.9412

cat("Cosine similarity between BCL6nB and CD4T:", round(cosine_sim, 4), "\n")

## Cosine similarity between BCL6nB and CD4T: 0.3414

cat("Pearson correlation between BCL6nB and CD4T:", round(pearson_cor, 4), "\n")

## Pearson correlation between BCL6nB and CD4T: 0.2188

# Create comparison table
correlation_comparison <- data.frame(
  Method = c("Graph Cross-Correlation (GCC)", "Cosine Similarity", "Pearson Correlation"),
  Value = c(round(gcc_result, 4), round(cosine_sim, 4), round(pearson_cor, 4)),
  Description = c("Frequency-domain correlation using graph structure",
                  "Geometric similarity measure",
                  "Linear correlation coefficient")
)

kable(correlation_comparison, 
      caption = "Cross-Cell Type Correlation Analysis: BCL6nB vs CD4T")

Cross-Cell Type Correlation Analysis: BCL6nB vs CD4T

Method	Value	Description
Graph Cross-Correlation (GCC)	0.9412	Frequency-domain correlation using graph structure
Cosine Similarity	0.3414	Geometric similarity measure
Pearson Correlation	0.2188	Linear correlation coefficient

BioGSP_Tutorial.R

Cross-ROI Analysis

cat("=== Demonstrating Cross-ROI Analysis ===\n")

## === Demonstrating Cross-ROI Analysis ===

# Compare BCL6- B Cell patterns across different ROIs
roi_comparison <- list()
for (roi in c("ROI_1", "ROI_2", "ROI_3")) {
  roi_data <- codex_toy_data %>% filter(ROI_num == roi)
  roi_bcl6nb <- roi_data %>% filter(Annotation5 == "BCL6- B Cell")
  
  cat("ROI", roi, "- BCL6- B Cell:", nrow(roi_bcl6nb), "\n")
  cat("  X range:", round(range(roi_bcl6nb$X_cent), 2), "\n")
  cat("  Y range:", round(range(roi_bcl6nb$Y_cent), 2), "\n")
  
  # Quick binning for comparison
  if (nrow(roi_bcl6nb) > 20) {  # Only if sufficient cells
    binned <- create_binned_data(roi_data, "BCL6- B Cell", x_bins = 15, y_bins = 15)
    sgwt_data <- prepare_sgwt_data(binned)
    
    # Quick SGWT analysis
    result <- SGWT(data.in = sgwt_data,
                   signal = "signal",
                   k = min(8, nrow(sgwt_data) - 1),
                   J = 2,  # Fewer scales for quick analysis
                   kernel_type = "mexican_hat",
                   k_fold = 2,
                   return_all = TRUE)
    
    roi_comparison[[roi]] <- list(
      cells = nrow(roi_bcl6nb),
      rmse = result$reconstruction_error,
      active_bins = sum(binned$cell_count > 0)
    )
    
    cat("  SGWT RMSE:", round(result$reconstruction_error, 4), "\n")
    cat("  Active bins:", sum(binned$cell_count > 0), "out of 225\n\n")
  }
}

## ROI ROI_1 - BCL6- B Cell: 189 
##   X range: 28.14 59.85 
##   Y range: 76.06 106.72 
## Building graph from spatial coordinates...
## Computing Laplacian and eigendecomposition...
## Auto-generated scales: 0.5156, 0.2578 
## Performing SGWT decomposition...
## Reconstruction RMSE: 0.1359 
##   SGWT RMSE: 0.1359 
##   Active bins: 109 out of 225
## 
## ROI ROI_2 - BCL6- B Cell: 231 
##   X range: 49.92 87.46 
##   Y range: 13.05 46.97 
## Building graph from spatial coordinates...
## Computing Laplacian and eigendecomposition...
## Auto-generated scales: 0.5156, 0.2578 
## Performing SGWT decomposition...
## Reconstruction RMSE: 0.099712 
##   SGWT RMSE: 0.0997 
##   Active bins: 99 out of 225
## 
## ROI ROI_3 - BCL6- B Cell: 284 
##   X range: 36.21 95.79 
##   Y range: 16.7 70.47 
## Building graph from spatial coordinates...
## Computing Laplacian and eigendecomposition...
## Auto-generated scales: 0.5156, 0.2578 
## Performing SGWT decomposition...
## Reconstruction RMSE: 0.074979 
##   SGWT RMSE: 0.075 
##   Active bins: 96 out of 225

# Create comparison table
if (length(roi_comparison) > 0) {
  comparison_df <- data.frame(
    ROI = names(roi_comparison),
    BCL6_B_Cells = sapply(roi_comparison, function(x) x$cells),
    Active_Bins = sapply(roi_comparison, function(x) x$active_bins),
    SGWT_RMSE = round(sapply(roi_comparison, function(x) x$rmse), 4)
  )
  
  kable(comparison_df, caption = "Cross-ROI BCL6- B Cell Analysis Comparison")
}

Cross-ROI BCL6- B Cell Analysis Comparison

	ROI	BCL6_B_Cells	Active_Bins	SGWT_RMSE
ROI_1	ROI_1	189	109	0.1359
ROI_2	ROI_2	231	99	0.0997
ROI_3	ROI_3	284	96	0.0750

cat("Cross-ROI analysis demonstrates how spatial patterns vary across regions!\n")

## Cross-ROI analysis demonstrates how spatial patterns vary across regions!

BioGSP_Tutorial.R

Summary

cat("=== SGWT Analysis Summary ===\n\n")

## === SGWT Analysis Summary ===

cat("Dataset: Multi-ROI Toy CODEX spatial data\n")

## Dataset: Multi-ROI Toy CODEX spatial data

cat("- Total cells analyzed:", nrow(codex_toy_data), "\n")

## - Total cells analyzed: 18604

cat("- Number of ROIs:", length(unique(codex_toy_data$ROI_num)), "\n")

## - Number of ROIs: 16

cat("- ROI_1 focused analysis (30x30 grid):", 900, "bins\n")

## - ROI_1 focused analysis (30x30 grid): 900 bins

cat("- BCL6nB active bins (ROI_1):", sum(bcl6nb_binned$cell_count > 0), "\n")

## - BCL6nB active bins (ROI_1): 157

cat("- CD4T active bins (ROI_1):", sum(cd4t_binned$cell_count > 0), "\n")

## - CD4T active bins (ROI_1): 171

cat("\nSGWT Analysis Results:\n")

## 
## SGWT Analysis Results:

cat("- BCL6nB reconstruction RMSE:", round(sgwt_bcl6nb$reconstruction_error, 6), "\n")

## - BCL6nB reconstruction RMSE: 0.11396

cat("- CD4T reconstruction RMSE:", round(sgwt_cd4t$reconstruction_error, 6), "\n")

## - CD4T reconstruction RMSE: 0.149172

cat("- Number of scales:", length(sgwt_bcl6nb$parameters$scales), "\n")

## - Number of scales: 3

cat("- Kernel type:", sgwt_bcl6nb$parameters$kernel_type, "\n")

## - Kernel type: mexican_hat

cat("\nCross-Correlation Results:\n")

## 
## Cross-Correlation Results:

cat("- Graph Cross-Correlation:", round(gcc_result, 4), "\n")

## - Graph Cross-Correlation: 0.9412

cat("- Cosine Similarity:", round(cosine_sim, 4), "\n")

## - Cosine Similarity: 0.3414

cat("- Pearson Correlation:", round(pearson_cor, 4), "\n")

## - Pearson Correlation: 0.2188

if (abs(gcc_result) > 0.3) {
  cat("\n✓ Moderate to strong spatial correlation detected between cell types\n")
} else {
  cat("\n✓ Weak spatial correlation between cell types\n")
}

## 
## ✓ Moderate to strong spatial correlation detected between cell types

cat("\n✓ Tutorial completed successfully using the BioGSP package!\n")

## 
## ✓ Tutorial completed successfully using the BioGSP package!

BioGSP_Tutorial.R

Conclusion

This tutorial demonstrated the key capabilities of the BioGSP package:

1. Data Integration: Using built-in toy CODEX data for reproducible analysis
2. Kernel Families: Comparing Mexican hat, Meyer, and heat kernel families
   
3. Spatial Analysis: Multi-scale decomposition of spatial cell patterns
4. Graph Analysis: Construction and validation of spatial graphs
5. Correlation Analysis: Advanced measures for cross-cell type relationships

The BioGSP package provides a comprehensive framework for spatial analysis of multiplexed imaging data, enabling researchers to uncover multi-scale spatial patterns and relationships in complex tissue environments.

For more information, see ?codex_toy_data for dataset details and explore other BioGSP functions using help(package = "BioGSP").